1. Field of the Invention
The present invention is a method of obtaining an anisotropic velocity model for proper depth imaging of seismic data.
2. Related Prior Art
The search for subsurface hydrocarbon deposits typically involves a multifaceted sequence of data acquisition, analysis, and interpretation procedures. The data acquisition phase involves use of an energy source to generate signals which propagate into the earth and reflect from various subsurface geologic structures. The reflected signals are recorded by a multitude of receivers on or near the surface of the earth, or in an overlying body of water. The received signals, which are often referred to as seismic traces, consist of amplitudes of acoustic energy which vary as a function of time, receiver position, and source position and, most importantly, vary as a function of the physical properties of the structures from which the signals reflect. The data analyst uses these traces along with a geophysical model to develop an image of the subsurface geologic structures.
Common Mid Point (CMP) stacking, also sometimes referred to as Common Depth Point or Common Reflection Point, CDP or CRP respectively) of seismic field data is well known. See for example U.S. Pat. No. 3,217,828 to Mendenhall et al, and U.S. Pat. No. 2,732,906 to Mayne, which are incorporated herein by reference as a teaching of the CMP technique. The U.S. Pat. No. 3,217,828 teaches two-dimensional (2-D) data processing where dipping earth layers are projected into a two-dimensional plane, perpendicular to the surface of the earth, along a designated line of profile.
Seismic data migration typically uses diffraction traveltimes from subsurface imaging points to the source and receiver locations to produce an image of the subsurface reflectors. The diffraction traveltimes are the seismic signal propagation times along raypaths from each imaging point to the source and receiver locations. The propagation times, which are usually plotted as diffraction traveltime curves, are used after appropriate preprocessing of the raw seismic data to generate an estimate of the correct location of the reflector.
The basic concepts of migration are discussed next with reference to FIG. 1. A portion of the surface of the earth is shown at 10. A sound source S may be offset from a receiver R by an offset distance 2 k with a midpoint at m0. A wavefield generated at S, may travel along ray path 12 to a reflecting point 14 on reflector 16 and reflected back along ray path 18 to receiver R. The amplitude of the wave field as a function of time, f(t), may be recorded on a time-scale recording, hereinafter referred to as trace. Seismic data are generally quantized as digital samples. In the processes next to be discussed, each sample of a seismic trace may be operated upon individually. In the interest of brevity, use of the collective term “seismic trace” in conjunction with a process means that every data sample of that trace has been individually subjected to the named process.
Because of the offset distance 2 k, the arrival time t of a particular wavelet along path S-14-R is greater than the travel time t0 of a wavelet that might have traveled along a direct path m0-14-m0. The time difference is termed normal moveout (NMO). The quantity f(t) may be corrected for NMO by the relation                               t          0          2                =                              t            k            2                    -                      (                          4              ⁢                                                k                  2                                /                                  V                  2                                                      )                                              (        1        )            where V represents the NMO velocity of the medium through which the wavelet traveled. The above relation is a small-offset approximation. It should be observed that for zero offset k and constant velocity V the reflection points for f(t0) such as 14 lie along a semicircle such as 19 centered about m0. Wavelet amplitude varies according to the inverse square law due to geometric spreading.
The hyperbolic time-distance relationship given by eq. (1) is illustrated in FIG. 2 where a plurality of traces 21a, 21b, . . . 21n are shown with a reflection event depicted by the curve 25. The general process of migration involves moving portions of seismic data along a diffraction curve such as 25 to an output location and summing the contributions from a large portion of the seismic section. The reflector in the subsurface may be viewed as being made up of a plurality of diffraction points. At the correct position of a diffractor, contributions from the plurality of migrated traces will sum coherently whereas at other locations, the summation will be incoherent and have close to zero amplitude; in this fashion an image of the subsurface can be reconstructed.
The migration process may be implemented in the time domain by the summation process or it may be implemented in the frequency domain. Such methods would be known to those versed in the art.
There is a vast array of literature indicating that the velocity of propagation of seismic waves in the subsurface is anisotropic, i.e., it depends upon the direction of propagation. There are numerous causes of such anisotropy. Backus showed that layering of isotropic materials having different shear modulii and with a layer thickness much shorter than the wavelength of seismic waves produces a transverse isotropy with a vertical symmetry axis (referred to hereafter VTI). In addition, shales exhibit transverse isotropy in their elastic properties even at very short seismic wavelengths. This has been attributed to intrinsic anisotropy. The cause of VTI is not pertinent to the present invention in that the present invention is a method for properly migrating seismic data in the presence of VTI regardless of the origin of the VTI.
A TI medium is characterized by five elastic modulii. These may be denoted by the tensor                     [                                                            c                11                                                                                      c                  11                                -                                  2                  ⁢                                      c                    66                                                                                                      c                13                                                    0                                      0                                      0                                                                                            c                  11                                -                                  2                  ⁢                                      c                    66                                                                                                      c                11                                                                    c                13                                                    0                                      0                                      0                                                                          c                13                                                                    c                13                                                                    c                33                                                    0                                      0                                      0                                                          0                                      0                                      0                                                      c                44                                                    0                                      0                                                          0                                      0                                      0                                      0                                                      c                44                                                    0                                                          0                                      0                                      0                                      0                                      0                                                      c                66                                                    ]                            (        2        )            where the modulus c11 defines the velocity of a horizontally propagating P-wave, c33 defines the velocity of a vertically propagating P-wave, c44 defines the velocity of a vertically propagating shear wave, and c66 defines the velocity of a horizontally propagating S-wave (shear wave) with horizontal polarization. These four parameters are determinable by making suitable measurements of P- and S-waves parallel to and perpendicular to the symmetry axis. The velocity is given by the square root of the ratio of the elastic modulus to the density. However, the parameter c13 is not readily determinable by measurements parallel to and perpendicular to the symmetry axis. Determination of c13 requires a measurement off the symmetry axis. It also turns out that this parameter c13 is of particular importance in reflection seismic prospecting.
Thomsen defined several parameters characterizing a VTI medium. One of these is the quantity ε defined by                     ɛ        =                                            c              11                        -                          c              33                                            2            ⁢                          c              33                                                          (        3        )            This parameter is directly related to the difference between the horizontal and vertical P-wave velocities. Another parameter defined by Thomsen is the quantity δ defined as                     δ        =                                                            (                                                      c                    13                                    +                                      c                    44                                                  )                            2                        -                                          (                                                      c                    33                                    -                                      c                    44                                                  )                            2                                            2            ⁢                                          c                33                            ⁡                              (                                                      c                    33                                    -                                      c                    44                                                  )                                                                        (        4        )            Based on the results derived by Thomsen, for a horizontal reflector in a TI medium, eq.(1) is valid for small offsets withtn=z/Vz  (5)andVnmo=V=Vz√{square root over (1+2δ)}  (6)where z is the depth of the reflector, Vz is the vertical velocity. The normal moveout velocity is thus dependent upon the parameter δ and the vertical velocity.
Alkhalifah & Tsvankin showed that the P-wave reflection moveout at large offsets is largely controlled by just two parameters. One is the Vnmo given by eq. (6) and the second is a parameter η defined by                     η        =                                            1              2                        ⁢                          (                                                                    V                    hor                    2                                                        V                    nmo                    2                                                  -                1                            )                                =                                    ɛ              -              δ                                      1              +                              2                ⁢                δ                                                                        (        7        )            whereVhor=Vnmo√{square root over (1+2η)}  (8)with the moveout equation given by                                                                                           t                  2                                ⁡                                  (                  x                  )                                            =                            ⁢                                                t                  o                  2                                +                                                      x                    2                                                        V                    nmo                    2                                                  -                                                      2                    ⁢                                                                                   ⁢                    η                    ⁢                                                                                   ⁢                                          x                      4                                                                                                  V                      nmo                      2                                        [                                                                                            t                          0                          2                                                ⁢                                                  V                          nmo                          2                                                                    +                                                                        (                                                      1                            +                                                          2                              ⁢                                                                                                                           ⁢                              η                                                                                )                                                ⁢                                                  x                          2                                                                                      ]                                                                                                                          =                            ⁢                                                t                  0                  2                                +                                                      x                    2                                                        V                    nmo                    2                                                  -                                                                            (                                                                        V                          hor                          2                                                -                                                  V                          nmo                          2                                                                    )                                        ⁢                                          x                      4                                                                                                  V                      nmo                      2                                        (                                                                                            t                          2                                                ⁢                                                  V                          nmo                          4                                                                    +                                                                        V                          ho                          2                                                ⁢                                                  x                          2                                                                                      )                                                                                                          (        9        )            and x=2 k in eq. (1) and FIG. 1.
As would be known to those versed in the art, conventional velocity analysis of seismic reflection data is based on use of semblance scans along hyperbolic moveout curves. Gretchka & Tsvankin disclose a method for inversion of seismic data using semblance scans along non-hyperbolic moveout curves using a modified version of eq. (9). The modification is based on empirical observations for simulated results wherein an offset range of x/z≅2 was used. The modified equation used by Gretchka & Tsvankin is                                           t            2                    ⁡                      (            x            )                          =                              t            0            2                    +                                    x              2                                      V              nmo              2                                -                                                    (                                                      V                    hor                    2                                    -                                      V                    nmo                    2                                                  )                            ⁢                              x                4                                                                    V                nmo                2                            (                                                                    t                    0                    2                                    ⁢                                      V                    nmo                    4                                                  +                                                      CV                    hor                    2                                    ⁢                                      x                    2                                                              )                                                          (        10        )            the modification being the empirical constant C in the denominator of the last term on the right hand side of eq. (10).
Accuracy of the subsurface image is of great importance in seismic prospecting: interpretation of the structure and lithology of the subsurface is preferably done on depth migrated sections. The accuracy of the image is closely related to the velocity model. Gretchka & Tsvankin do not address the interaction between iterative changes in the velocity model and subsequent changes in the accuracy of the depth image.
It is therefore desirable to have a method for determination and iterative updating of subsurface velocities in anisotropic earth formations. Such a method should preferably be efficient in terms of computer resources and should also preferably provide more accurate images of the earth formations. The present invention satisfies this need.